#ifndef _FSIP7_H_
#define _FSIP7_H_

#include "Equations/TwoPhaseEquation.hpp"

namespace Tuna {

  template <typename T, int Dim> class FSIP7;

  template <typename T>
  struct Typeinfo <FSIP7<T, 1> > {
    typedef T prec_t;
    enum { Dim = 1 };
  };

  template <typename T>
  struct Typeinfo <FSIP7<T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };

  template <typename T>
  struct Typeinfo <FSIP7<T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };

  
  template<typename Tprec, int Dim>
  class FSIP7 : public TwoPhaseEquation<FSIP7<Tprec, Dim> >
  {

    typedef TwoPhaseEquation<FSIP7<Tprec, Dim> > TP_FSIP7;
    
    using GeneralEquation< TP_FSIP7 >::aE;
    using GeneralEquation< TP_FSIP7 >::aW;
    using GeneralEquation< TP_FSIP7 >::aN;
    using GeneralEquation< TP_FSIP7 >::aS;
    using GeneralEquation< TP_FSIP7 >::aF;
    using GeneralEquation< TP_FSIP7 >::aB;
    using GeneralEquation< TP_FSIP7 >::aP;
    using GeneralEquation< TP_FSIP7 >::sp;
    using GeneralEquation< TP_FSIP7 >::dx;
    using GeneralEquation< TP_FSIP7 >::dy;
    using GeneralEquation< TP_FSIP7 >::dz;
    using GeneralEquation< TP_FSIP7 >::bi;
    using GeneralEquation< TP_FSIP7 >::ei;
    using GeneralEquation< TP_FSIP7 >::bj;
    using GeneralEquation< TP_FSIP7 >::ej;
    using GeneralEquation< TP_FSIP7 >::bk;
    using GeneralEquation< TP_FSIP7 >::ek;
    using GeneralEquation< TP_FSIP7 >::applyBoundaryConditions1D;
    using GeneralEquation< TP_FSIP7 >::applyBoundaryConditions2D;
    using GeneralEquation< TP_FSIP7 >::applyBoundaryConditions3D;

    using TP_FSIP7::S;
    using TP_FSIP7::phi_0;
    using TP_FSIP7::Srw;
    using TP_FSIP7::Sro;
    using TP_FSIP7::mu_w;
    using TP_FSIP7::mu_o;
    using TP_FSIP7::k;
    using TP_FSIP7::rank;		//nuevo
    using TP_FSIP7::injection;  //nuevo

  public:
    typedef Tprec prec_t;
    typedef typename TunaArray<prec_t, Dim >::huge ScalarField;
    
    FSIP7() : TwoPhaseEquation<FSIP7<prec_t, Dim > >() { }    
    ~FSIP7() { };
    
    inline bool calcCoefficients1D() { }; 
    inline bool calcCoefficients2D() { };
    inline bool calcCoefficients3D();
    inline void printInfo() { std::cout << " FSIP7 "; }
  };

  
/*
 *  Lineal for realtive permeability , Upwind for Sw
 *
template<typename Tprec, int Dim>
inline bool FSIP7<Tprec, Dim>::calcCoefficients1D () 
{
    static prec_t Sw_e, Sw_w;

    // Lineal 
    static prec_t mult_o = k / ( (1 - Srw - Sro) * mu_o * dx ) ;
    static prec_t mult_w = k / ( (1 - Srw - Sro) * mu_w * dx ) ;

    aE = 0.0; aW = 0.0; aP = 0.0; sp = 0.0;

    for (int i =  bi; i <= ei; ++i) {  
	  
	  // UpwindQ
      if ( phi_0(i+1) >= phi_0(i) )
	if      ( i == ei   ) Sw_e = S(i+1);
	else if ( i == ei-1 ) Sw_e = S(i+1) * 2.0 - S(i+2);
	else Sw_e = S(i+1) + 0.125 * ( 3 * S(i) - 2 * S(i+1) - S(i+2) );
      else
	if ( i == bi ) Sw_e = S(i);
	else Sw_e = S(i) + 0.125 * (3 * S(i+1) - 2 * S(i) - S(i-1) );

      if ( phi_0(i-1) >= phi_0(i) ) 
	if      ( i == bi   ) Sw_w = S(i-1);
	else if ( i == bi+1 ) Sw_w = S(i-1);
	else Sw_w = S(i-1) + 0.125 * (3 * S(i) - 2 * S(i-1) - S(i-2) );
      else
	if ( i == ei ) Sw_w = S(i);
	else Sw_w = S(i) + 0.125 * (3 * S(i-1) - 2 * S(i) - S(i+1) );
    
      
      // Lineal
      aE (i) = (1 - Sro - Sw_e) * mult_o + (Sw_e - Srw) * mult_w ;
      aW (i) = (1 - Sro - Sw_w) * mult_o + (Sw_w - Srw) * mult_w ;
      aP (i) = aE (i) + aW (i);      
    }
    applyBoundaryConditions1D();
    return 0;
}

/*
 *  Lineal for relative permeability , Upwind for Sw
 *
template<typename Tprec, int Dim>
inline bool FSIP7<Tprec, Dim>::calcCoefficients2D () 
{
  static prec_t Sw_e, Sw_w, Sw_n, Sw_s;

    // Lineal 
    static prec_t mult_o = k / ( (1 - Srw - Sro) * mu_o ) ;
    static prec_t mult_w = k / ( (1 - Srw - Sro) * mu_w ) ;
    static prec_t dx_dy = dx / dy;
    static prec_t dy_dx = dy / dx;
  
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; sp = 0.0;


    for (int i =  bi; i <= ei; ++i)
      for (int j = bj; j <= ej; ++j) {  
      // UpwindE
// -- X
	if ( phi_0(i+1, j) >= phi_0(i, j) )
	  if      ( i == ei   ) Sw_e = S(i+1, j);
	  else if ( i == ei-1 ) Sw_e = S(i+1, j) * 2.0 - S(i+2, j);
	  else                  Sw_e = S(i+1, j) + 0.125 * ( 3 * S(i,j) - 2 * S(i+1, j) - S(i+2, j) );
      else
	if ( i == bi )
	  Sw_e = S(i, j);
	else
	  Sw_e = S(i, j) + 0.125 * (3 * S(i+1, j) - 2 * S(i, j) - S(i-1, j) );
	
	if ( phi_0(i-1, j) >= phi_0(i, j) ) 
	  if      ( i == bi   ) Sw_w = S(i-1, j);
	  else if ( i == bi+1 ) Sw_w = S(i-1, j) ;
	  else                  S(i-1, j) + 0.125 * (3 * S(i,j) - 2 * S(i-1, j) - S(i-2, j) );
      else
	if ( i == ei ) Sw_w = S(i, j);
	else
	  Sw_w = S(i,j) + 0.125 * (3 * S(i-1, j) - 2 * S(i,j) - S(i+1, j) ); 

// -- Y
      if ( phi_0(i, j+1) >= phi_0(i,j) )
	if      ( j == ej  ) Sw_n = S(i, j+1);
	else if ( j == ej-1 ) Sw_n = S(i, j+1) * 2.0 - S(i, j+2);
	else Sw_n = S(i, j+1) + 0.125 * ( 3 * S(i, j) - 2 * S(i, j+1) - S(i, j+2) );
      else
	if ( j == bj ) 
		Sw_n = S(i,j);
	else Sw_n = S(i,j) + 0.125 * (3 * S(i, j+1) - 2 * S(i, j) - S(i, j-1) );

      if ( phi_0(i, j-1) >= phi_0(i, j) ) 
	if      ( j == bj   ) Sw_s = S(i, j-1);
	else if ( j == bj+1 ) Sw_s = S(i, j-1);
	else Sw_s = S(i, j-1) + 0.125 * (3 * S(i,j) - 2 * S(i, j-1) - S(i, j-2) );
      else
	if ( j == ej ) Sw_s = S(i, j);
	else Sw_s = S(i, j) + 0.125 * (3 * S(i, j-1) - 2 * S(i, j) - S(i, j+1) );
      
      // Lineal
	aE (i, j) = ( (1 - Sro - Sw_e) * mult_o +
		      (Sw_e - Srw) * mult_w ) * dy_dx;
	aW (i, j) = ( (1 - Sro - Sw_w) * mult_o + 
		      (Sw_w - Srw) * mult_w ) * dy_dx;
	aN (i, j) = ( (1 - Sro - Sw_n) * mult_o + 
		      (Sw_n - Srw) * mult_w ) * dx_dy;
	aS (i, j) = ( (1 - Sro - Sw_s) * mult_o + 
		      (Sw_s - Srw) * mult_w ) * dx_dy;
	aP (i, j) = aE (i, j) + aW (i, j) + aN (i, j) + aS (i, j);
      }
    applyBoundaryConditions2D();
    
    //nuevo
        // Source for injection : 
    // (\frac{\bar{q}_w}{\rho_w} + \frac{\bar{q}_w}{\rho_w}) * \Delta V
    sp (bi, bj) +=  injection  ;

/*
    sp (ei, ej) +=  injection  ;
    sp (bi, ej) +=  injection  ;
    sp (ei, bj) +=  injection  ;
*
    //Sink for extraction :
    sp (ei, ej) -= injection  ;

//    sp (ei/2+1, ej/2+1) -= 4 * injection  ; 
//hasta aca nuevo
    
    return 0;
}

/*
 *  Lineal for relative permeability , Upwind for Sw
 */
template<typename Tprec, int Dim>
inline bool FSIP7<Tprec, Dim>::calcCoefficients3D () 
{
    static prec_t Sw_e, Sw_w, Sw_n, Sw_s, Sw_f, Sw_b;

    // Lineal 
    static prec_t mult_o = k / ( (1 - Srw - Sro) * mu_o ) ;
    static prec_t mult_w = k / ( (1 - Srw - Sro) * mu_w ) ;
    static prec_t dxdz_dy = dx * dz / dy;
    static prec_t dydz_dx = dy * dz / dx;
    static prec_t dxdy_dz = dx * dy / dz;

    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aF = 0.0; aB = 0.0;
    aP = 0.0; sp = 0.0;

    for (int ki = bk; ki <= ek; ++ki) 
      for (int i =  bi; i <= ei; ++i) 
	for (int j = bj; j <= ej; ++j) {
      
      // UpwindQ
// -- X
	  if ( phi_0(i+1, j, ki) >= phi_0(i, j, ki) )
	    if      ( i == ei   ) Sw_e = S(i+1, j, ki);
	    else if ( i == ei-1 ) Sw_e = S(i+1, j, ki) * 2.0 - S(i+2, j, ki);
	    else Sw_e = S(i+1, j, ki) + 0.125 * ( 3 * S(i, j, ki) - 2 * S(i+1, j, ki) - S(i+2, j, ki) );
	  else
	    if ( i == bi ) Sw_e = S(i, j, ki);
	    else Sw_e = S(i, j, ki) + 0.125 * (3 * S(i+1, j, ki) - 2 * S(i, j, ki) - S(i-1, j, ki) );

	  
	  if ( phi_0(i-1, j, ki) >= phi_0(i, j, ki) ) 
	    if      ( i == bi   ) Sw_w = S(i-1, j, ki);
	    else if ( i == bi+1 ) Sw_w = S(i-1, j, ki);
	    else Sw_w = S(i-1, j, ki) + 0.125 * (3 * S(i, j, ki) - 2 * S(i-1, j, ki) - S(i-2, j, ki) );
	  else
	    if ( i == ei ) Sw_w = S(i, j, ki);
	    else Sw_w = S(i, j, ki) + 0.125 * (3 * S(i-1, j, ki) - 2 * S(i, j, ki) - S(i+1, j, ki) );

// -- Y
	  if ( phi_0(i, j+1, ki) >= phi_0(i, j, ki) )
	    if      ( j == ej   ) Sw_n = S(i, j+1, ki);
	    else if ( j == ej-1 ) Sw_n = S(i, j+1, ki) * 2.0 - S(i, j+2, ki);
	    else Sw_n = S(i, j+1, ki) + 0.125 * ( 3 * S(i, j, ki) - 2 * S(i, j+1, ki) - S(i, j+2, ki) );
	  else
	    if ( j == bj ) Sw_n = S(i, j, ki);
	    else Sw_n = S(i, j, ki) + 0.125 * (3 * S(i, j+1, ki) - 2 * S(i, j, ki) - S(i, j-1, ki) );

	  
	  if ( phi_0(i, j-1, ki) >= phi_0(i, j, ki) ) 
	    if      ( j == bj   ) Sw_s = S(i, j-1, ki);
	    else if ( j == bj+1 ) Sw_s = S(i, j-1, ki);
	    else Sw_s = S(i, j-1, ki) + 0.125 * (3 * S(i, j, ki) - 2 * S(i, j-1, ki) - S(i, j-2, ki) );
	  else
	    if ( j == ej ) Sw_s = S(i, j, ki);
	    else Sw_s = S(i, j, ki) + 0.125 * (3 * S(i, j-1, ki) - 2 * S(i, j, ki) - S(i, j+1, ki) );

// -- Z
	  if ( phi_0(i, j, ki+1) >= phi_0(i, j, ki) )
	    if      ( ki == ek   ) Sw_f = S(i, j, ki+1);
	    else if ( ki == ek-1 ) Sw_f = S(i, j, ki+1) * 2.0 - S(i, j, ki+2);
	    else Sw_f = S(i, j, ki+1) + 0.125 * ( 3 * S(i, j, ki) - 2 * S(i, j, ki+1) - S(i, j, ki+2) );
	  else
	    if ( ki == bk ) Sw_f = S(i, j, ki);
	    else Sw_f = S(i, j, ki) + 0.125 * (3 * S(i, j, ki+1) - 2 * S(i, j, ki) - S(i, j, ki-1) );

	  
	  if ( phi_0(i, j, ki-1) >= phi_0(i, j, ki) ) 
	    if      ( ki == bk   ) Sw_f = S(i, j, ki-1);
	    else if ( ki == bk+1 ) Sw_f = S(i, j, ki-1);
	    else Sw_s = S(i, j, ki-1) + 0.125 * (3 * S(i, j, ki) - 2 * S(i, j, ki-1) - S(i, j, ki-2) );
	  else
	    if ( k == ek ) Sw_b = S(i, j, ki);
	    else Sw_b = S(i, j, ki) + 0.125 * (3 * S(i, j, ki-1) - 2 * S(i, j, ki) - S(i, j, ki+1) );
	  
	  // Lineal
	  aE (i, j, ki) = ( (1 - Sro - Sw_e) * mult_o + 
			   (Sw_e - Srw) * mult_w ) * dydz_dx;
	  aW (i, j, ki) = ( (1 - Sro - Sw_w) * mult_o + 
			   (Sw_w - Srw) * mult_w ) * dydz_dx;
	  aN (i, j, ki) = ( (1 - Sro - Sw_n) * mult_o + 
			   (Sw_n - Srw) * mult_w ) * dxdz_dy;
	  aS (i, j, ki) = ( (1 - Sro - Sw_s) * mult_o + 
			   (Sw_s - Srw) * mult_w ) * dxdz_dy;
	  aF (i, j, ki) = ( (1 - Sro - Sw_f) * mult_o + 
			   (Sw_f - Srw) * mult_w ) * dxdy_dz;
	  aB (i, j, ki) = ( (1 - Sro - Sw_b) * mult_o + 
			   (Sw_b - Srw) * mult_w ) * dxdy_dz;	  
	  aP (i, j, ki) = aE (i, j, ki) + aW (i, j, ki) + 
	                  aN (i, j, ki) + aS (i, j, ki) + 
	                  aF (i, j, ki) + aB (i, j, ki);
	}        
    applyBoundaryConditions3D();

    // Source for injection : 
    // (\frac{\bar{q}_w}{\rho_w} + \frac{\bar{q}_w}{\rho_w}) * \Delta V
    sp (bi, bj, bk) +=  injection  ;

/*
    sp (ei, ej) +=  injection  ;
    sp (bi, ej) +=  injection  ;
    sp (ei, bj) +=  injection  ;
*/
    //Sink for extraction :
    sp (ei, ej, ek) -= injection  ;

//    sp (ei/2+1, ej/2+1) -= 4 * injection  ;

    return 0;
}

} // Tuna namespace


#endif //_FSIP7_H_
